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Solidification 



O 

o 
o 

(N 
X> 

PL, 



O 

CO 



Mathis Plapp and Alain Karma 
Physics Department and Center for Interdisciplinary Research on Complex Systems, 
Northeastern University, Boston, Massachusetts 02115 
(January 11, 2000) 

We present a novel hybrid computational method to simulate accurately dendritic solidification in 
the low undercooling limit where the dendrite tip radius is one or more orders of magnitude smaller 
than the characteristic spatial scale of variation of the surrounding thermal or solutal diffusion field. 
The first key feature of this method is an efficient multiscale diffusion Monte-Carlo (DMC) algorithm 
which allows off-lattice random walkers to take longer and concomitantly rarer steps with increasing 
distance away from the solid- liquid interface. As a result, the computational cost of evolving the 
large scale diffusion field becomes insignificant when compared to that of calculating the interface 
evolution. The second key feature is that random walks are only permitted outside of a thin liquid 
layer surrounding the interface. Inside this layer and in the solid, the diffusion equation is solved 
using a standard finite-difference algorithm that is interfaced with the DMC algorithm using the 
local conservation law for the diffusing quantity. Here we combine this algorithm with a previously 
developed phase-field formulation of the interface dynamics and demonstrate that it can accurately 
simulate three-dimensional dendritic growth in a previously unreachable range of low undercoolings 
that is of direct experimental relevance. 
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I. INTRODUCTION 

Diffusion-limited pattern formation, which leads to the 
spontaneous emergence of complex branched structures, 
occurs in numerous contexts. A few examples include 
dendritic solidification , electrochemical deposition 0] 
and corrosion, and the growth of bacterial colonies [p| . 
Two distinct length scales are typically involved in this 
class of problems: one that characterizes the pattern it- 
self, such as the thickness of a branch, and one that char- 
acterizes the diffusion field associated with the transport 
of heat or matter. In many cases, these two scales are 
vastly different. For example, in solidification, the decay 
length of the thermal or solutal field ahead of a growing 
dendrite (in a pure or alloy melt) can be one to three or- 
ders of magnitude larger than the tip radius of one of its 
primary branches. Non-trivial pattern formation dynam- 
ics can be expected to occur on all intermediate scales. 
This poses a serious challenge for numerical simulations 
since a precise integration of the equations of motion on 
the pattern scale requires a good resolution of the interfa- 
cial region, and such a resolution is completely inefficient 
(i.e. much too fine) to treat the large scale diffusion field. 
Therefore, in order to retain this precision on the small 
scale and, at the same time, simulate the pattern evo- 
lution on sufficiently large length and time scales, it is 
necessary to use some form of multiscale algorithm. 

Multi-grid and finite element methods with non- 
uniform meshing represent one possible solution for this 
type of problems. Their application, however, in the con- 
text of growth simulations faces the additional difficulty 
of a moving interface, which implies that the structure 
of the simulation grid has to be dynamically adapted. 
For the classic problem of dendritic crystal growth, sev- 



eral multi-grid Q or adaptive meshing algorithms || 
have been proposed in recent years. The most precise 
to date is the method of Provatas et al. which uses the 
phase-field model on a regular grid on the scale of the 
dendrite, whereas the diffusion field is integrated on an 
adaptive mesh using finite element techniques [[| . While 
this method appears to be promising, it has yet to be 
implemented in three dimensions where the difficulty of 
adaptive meshing becomes significantly enhanced. 

We present in this paper an alternative solution to 
solve this computational challenge and we illustrate its 
application in the context of the dendritic crystalliza- 
tion of a pure substance from its undercooled melt, even 
though this algorithm can be applied to any diffusion- 
limited growth problem for which an explicit solver of 
the interface dynamics is available. The idea is to use a 
hybrid approach. The interface dynamics is treated using 
deterministic equations of motion, in particular those of 
the phase-field model for the dendritic growth problem 
considered here. On the other hand, the large-scale diffu- 
sion field is represented by an ensemble of off-lattice ran- 
dom walkers and is evolved using a Diffusion Monte Carlo 
(DMC) algorithm. The two solutions are connected at 
some distance from the moving interface. The key point 
for rendering our method efficient is that we use random 
walkers which dynamically adapt the average length of 
their random steps. Far from the interface, the walkers 
can make large jumps and hence be updated only rarely 
without affecting the quality of the solution near the 
growing interface. In some sense, our method can be seen 
as an "adaptive grid algorithm without grid" . The DMC 
algorithm and the connection between deterministic and 
stochastic parts are rather simple and straightforward to 
implement in both two and three dimensions, both on 
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single-processor and parallel architectures. We demon- 
strate in this paper that our method is precise, robust, 
and reliable, and hence constitutes a powerful alternative 
to state of the art adaptive meshing techniques. Tech- 
nically, the algorithm bears many similarities to quan- 
tum Monte Carlo methods [Q . It is therefore remarkable 
that the gap between mesoscopic and macroscopic length 
scales can be bridged using a method borrowed from mi- 
croscopic physics in an interfacial pattern formation con- 
text, which was not a priori obvious to us at the start of 
this investigation. 

Our algorithm builds on ideas of earlier random walk 
algorithms for simulating pattern formation during vis- 
cous fingering ||,|| and solidification JToj-p^l, but intro- 
duces two essential new features. Firstly, random walks 
with variable step size have been used previously in sim- 
ulations of large-scale diffusion- limited aggregation jli} , 
but only one walker at a time was simulated, and hence 
the time variable did not explicitly appear in the treat- 
ment of the walkers. In the present diffusive case, the 
memory of the past history, stored in the diffusion field, is 
essential to the problem. Our DMC algorithm works with 
a whole ensemble of walkers in "physical" time and hence 
constitutes a true multiscale solver for the full diffusion 
problem. Secondly, the algorithms mentioned above use 
a lattice both to evolve the walkers and to represent the 
position of the interface by the bonds between occupied 
(solid) and empty (liquid) sites. Walkers are created or 
absorbed directly at this interface. The discretization of 
space and the stochastic creation and absorption of walk- 
ers make it difficult to control accurately the interfacial 
anisotropy and the noise that both play a crucial role 
in dendritic evolution [ pT|p^ ]. Consequently, the algo- 
rithms aimed at describing dendritic growth fl(i| , |ll] . |l3| . 
while correctly reproducing all the qualitative features of 
the growth process, are unable to yield quantitative re- 
sults that can be tested against experiments. We solve 
both problems by creating and absorbing walkers not at 
the solid-liquid interface, but at a "conversion boundary" 
at some fixed distance from the interface. This means 
that the stochastic representation of the diffusion and the 
motion of the interface can be treated separately, which 
allows us to evolve the interface accurately by the phase- 
field method using a finite difference representation of 
controlled precision. At the same time, the stochastic 
noise created by the DMC algorithm is rapidly damped 
by the deterministic diffusion in the "buffer layer" be- 
tween the conversion boundary and the solid-liquid in- 
terface, and hence the amplitude of the fluctuations at 
the solid-liquid interface can be reduced to a prescribed 
level without much cost in computation time by increas- 
ing the thickness of the buffer layer. This is an impor- 
tant issue for simulations of dendritic growth, because 
the amplification of microscopic fluctuations of the inter- 
face is believed to be the main cause for the formation 
of secondary dendrite branches , and it is well known 
that numerical noise can lead to the formation of spuri- 
ous sidebranches in simulations. Consequently, we have 



to demonstrate that the walker noise of our algorithm 
can be reduced to a level that does not affect the pattern 
evolution. 

Another benefit of the buffer layer is that it makes the 
algorithm very versatile. Away from the interface, only 
the standard diffusion equation has to be solved. There- 
fore, the DMC part of the algorithm and the conversion 
process between deterministic and stochastic solutions 
are completely independent from the method used for 
simulating the interface dynamics, and can easily be car- 
ried over to other free boundary problems. 

The purpose of the present paper is to describe the al- 
gorithm in detail and to demonstrate its reliability and 
precision by benchmark simulations. Some results con- 
cerning three-dimensional crystal growth at low under- 
coolings have already been presented elsewhere ]l7| , ^8t , 
and hence we will focus here on the computational as- 
pects of the problem. Section II contains a brief intro- 
duction to dendritic solidification and the basic equations 
of motion, and describes the phase-field method. In sec- 
tion III, the DMC algorithm and its interfacing with the 
phase-field equations are described in detail. In section 
IV, we present results of benchmark simulations, assess 
the efficiency of the code and the magnitude of numer- 
ical noise, and present simulations of three-dimensional 
dendritic growth. Section V contains a conclusion and 
the outline of future work. 



II. DENDRITIC GROWTH AND THE 
PHASE-FIELD METHOD 

When a crystal grows from an undercooled melt, it 
develops into an intricate branched structure, called a 
dendrite. This phenomenon has been of central impor- 
tance to the understanding of spontaneous pattern for- 
mation during phase transformations and the emergence 
of branched structures [f9|-|2l|l . In addition, it is of con- 
siderable practical interest, because dendrites form dur- 
ing the solidification of many commercially important al- 
loys and influence the mechanical properties of the fin- 
ished material. 

We will focus on the dendritic solidification of a pure 
substance from its homogeneously undercooled melt, 
starting from a single supercritical nucleus [^2| . This 
situation is well described by the symmetric model of so- 
lidification, which assumes that the diffusivity and ther- 
mophysical quantities such as the specific heat and the 
density are equal for the solid and the liquid phases. Dur- 
ing the growth of the crystal, the latent heat of melting 
is released, and in the absence of convection, the growth 
becomes limited by the diffusion of heat away from the 
growing dendrite. The state of the system at any time is 
described by the temperature field T(x, t) and the shape 
T(t) of the boundary between solid and liquid. It is cus- 
tomary to define a dimensionless temperature field 
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u(x, t) 



T(x, t) - T„ 
L/c p 



(1) 



where L and c p are the latent heat of melting and the 
specific heat, respectively, and T m is the melting temper- 
ature. In terms of this field, the equations of motion of 
the symmetric model are 



d t u = DV 2 u, 
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where D is the thermal diffusivity, v n is the normal ve- 
locity of the interface, and h is the unit normal vector of 
the surface T pointing towards the liquid. The diffusion 
equation, Eq. (^) is valid everywhere (in the liquid and 
in the solid) except on the surface T. The Stefan condi- 
tion, Eq. (|^), valid on T(t), expresses the conservation 
of enthalpy at the moving phase boundary. Here, Vm|s 
and Vu\l denote the limits of the temperature gradient 
when r is approached from the solid and the liquid side, 
respectively, and the equation states that the local heat 
flux at the interface must be equal to the latent heat 
generated or consumed during the phase transformation; 
v n is positive if the solid grows (i.e. freezing). The di- 
mensionless temperature at the interface ur is given by 
the generalized Gibbs-Thomson condition Eq. (|J). The 
first term on the right hand side (RHS) is the anisotropic 
form of the local equilibrium condition (Gibbs-Thomson 
condition) which relates the temperature to the curva- 
ture of the interface and the anisotropic surface tension 
7(n) = 7oa(n). For a crystal with cubic symmetry in 
three dimensions, the anisotropy function a(n) is usually 
written as 



a(h) = (1 — 3e4) 
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where £4 is the anisotropy parameter. Note that in two 
dimensions (d — 2), this expression reduces to 



a(6) = 1 + e 4 cos(46»), 



(6) 



where 9 is the angle between the normal and one of the 
axes of symmetry. On the RHS of Eq. (^) , 
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loT m c p 
L 2 
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is the capillary length, d is the spatial dimension, 8i 
are the angles between the normal n and the two local 
principal directions on T, and Ri are the principal radii 
of curvature. Finally, the second term on the RHS of 
Eq. (Eh describes the shift of the interface temperature 



due to molecular attachment kinetics, and (3(n) is the 
orientation-dependent linear kinetic coefficient. Kinetic 
effects are believed to be small for the range of solidifi- 
cation speeds of interest here. We will therefore focus on 
the case where the interface kinetics vanish (f3(n) = 0), 
which corresponds to local equilibrium at the interface. 
In this case, the physical length and time scales are set 
by the capillary length and the diffusivity, and the con- 
trol parameters of the problem are the anisotropy €4 and 
the dimensionless undercooling 



T m — Tq 
L/cp 



(8) 



where To is the initial temperature, T(x, 0) = To, which 
provides the thermodynamic driving force for solidifica- 
tion. We assume that the dendrite grows into an infinite 
volume of liquid, and hence u(x,t) — > — A as |x| — > 00 Vi. 
Typical experimental values for A range from 0.001 to 
0.1. The length scales involved in the problem are (i) 
the capillary length do , (ii) a typical scale of the pattern 
such as the radius of curvature at a tip p, and (iii) the 
length scale of the diffusion field Id- To fix the ideas, let 
us consider the measurements of Rubinstein and Glicks- 
man on pivalic acid (PVA) p3[ . For a dimensionless un- 
dercooling of A = 0.075, p = 8.5 fxm, and the speed of 
the tips is v = 390 fim/s, which gives a diffusion length 
Id = 2D/v = 0.38 mm, whereas do = 3.8 nm. The mul- 
tiscale character of this situation is obvious: Id and do 
differ by five orders of magnitude, and Id is forty times 
larger than p. These ratios become even larger for lower 
undercoolings. 

The above equations constitute a notoriously difficult 
free boundary problem. To simplify the task, theoreti- 
cal and numerical efforts first concentrated on the treat- 
ment of a single needle crystal growing at constant veloc- 
ity. This situation can be treated by boundary integral 
methods |2(J, which are exact in two dimensions (2-d) 
but have remained approximate in three dimensions (3- 
d). More recently, time-dependent methods have been 
developed to describe the full growth dynamics p9| . 
Of those, the phase-field method seems presently the 
most compact and precise approach. We use a recent effi- 
cient formulation of this method, which has been bench- 
marked against boundary integral calculations [^9). An 
"order parameter" , or phase-field ip{x, t) is introduced, 
which is an indicator field distinguishing the solid (-0 = 1) 
and the liquid (ip — —1) phase. The two-phase system is 
described by a free energy functional of Ginzburg-Landau 
type, 



dV[W 2 (n)\ViP\ 2 +/(</>, u)}, 



(9) 



where W{n) is the orientation-dependent interface thick- 
ness, i.e. the spatial scale on which the phase-field varies 
smoothly between its equilibrium values if) = ±1, and 
f(ip,u) is the free energy density. The equations of mo- 
tion are 
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(10) 



where SF/Sip denotes the functional derivative, and 



d t u = DW 2 u+-d t ^. 



(11) 



The phase-field relaxes to its local minimum free energy 
configuration, which depends on the local temperature 
field, with an orientation-dependent relaxation time r(n). 
The diffusion equation contains a source term to account 
for the latent heat released or consumed during the phase 
transformation. For a suitable choice of the functions 
/(i/>, u), W{n) and r(n), these equations reduce precisely 
to the free boundary problem given by Eqs. @ to (|J) in 
the limit where the interface thickness is small compared 
to the radii of curvature . A brief description of the 
model used for our simulations and its relation to the 
macroscopic free boundary problem is given in the ap- 
pendix. The key point is that the phase-field equations 
of motion are partial differential equations which can be 
integrated on a regular grid on the scale of the dendrite, 
without knowing explicitly where the solid-liquid inter- 
face is located. The phase field rapidly decays to its equi- 
librium values tp — ±1 away from the interface. There- 
fore, well within the bulk phases, Eq. (|l(]) becomes trivial 
and Eq. (11) reduces to the ordinary diffusion equation. 



III. DIFFUSION MONTE CARLO ALGORITHM 
A. Outline 

Our goal is to combine the precision of the phase-field 
method and the efficiency of a DMC treatment for the 
diffusion field. This is achieved by dividing the simula- 
tion domain into an "inner" and an "outer" region as 
shown in Fig. [l]. In the inner region, consisting of the 
growing structure and a thin "buffer layer" of liquid, we 
integrate the phase-field equations described above. In 
the outer region, the diffusion field is represented by an 
ensemble of random walkers. Walkers are created and ab- 
sorbed at the boundary between inner and outer domains 
at a rate which is proportional to the local diffusion flux. 
The value of the diffusion field in the outer domain is 
related to the local density of walkers, and the bound- 
ary conditions for the integration in the inner region are 
obtained by averaging this density over coarse-grained 
boxes close to the boundary. We will now describe in de- 
tail the DMC algorithm for the evolution of the random 
walkers and the connection of the two solutions. 

Let us start by recalling some well-known facts about 
random walkers. Consider first a single point particle 
performing a Brownian motion in continuous space and 
time. The conditional probability P(x' , t'\x, t) of finding 
the particle at position x 1 at time t' , given that it started 
from position x at time t, is identical to the diffusion 
kernel, 




FIG. 1. Simulation of two-dimensional dendritic growth for 
a dimensionless undercooling A = 0.1 and a surface ten- 
sion anisotropy £4 = 0.025. The solid line is the solid-liquid 
interface T, the dashed line is the conversion boundary V 1 
between the inner (deterministic) and outer (stochastic) do- 
mains, and the dots show the positions of random walkers 
(only one walker out of 50 is shown for clarity). 



P{x',t'\x, t) 



[A-nD{t> - t)] d/2 



exp 



— \x — XI 



W(t' - 1) 



(12) 



where D is the diffusion coefficient and d is the spatial 
dimension. This kernel satisfies the well-known convolu- 
tion relation 



P(x",t"\x,t)= J P(x",t"\x',t')P(x\t , \x,t)dx' 

Vt<t'<t". (13) 

Therefore, a realization of a random walk, i.e. the posi- 
tion of a walker as a function of time, represented by a 
time-dependent vector of real numbers x(t), can be ob- 
tained on a computer by successive steps. The position 
of the walker is updated following the scheme 



x(t + r) = x(t) + 



(14) 



where the components of the random vector £ are inde- 
pendent Gaussian random variables of unit variance. The 
time increment r (not to be confused with the phase-field 
relaxation time r(n) defined in the preceding section) and 
the step size I must satisfy the relation 



I 2 

— = 2D. 

T 



(15) 



Since time is continuous and Eq. ( |L3f ) is not restricted to 
t" — if = t' — t, successive steps may have different time 
increments (and concomitantly use different step lengths) 
if Eq. dT5l) is satisfied for each update. 



4 



The basic idea of Diffusion Monte Carlo simulations 
is to sample many realizations of diffusion paths. The 
density of random walkers then satisfies a stochastic dif- 
ferential equation which converges to the deterministic 
diffusion equation in the limit of an infinite number of 
walkers. A density of walkers can be defined by a suit- 
able coarse-graining procedure on a scale L cg , i.e. by 
dividing space into cells of volume L d cg and counting the 
number of walkers within each cell. If the coarse-graining 
length is chosen larger than the average step length £, this 
density evolves smoothly on the scale of L cg over times of 
order /D, the time for one walker to diffuse through 
a coarse cell. 

From the above considerations, it is clear that the char- 
acteristic length and time scales that can be resolved by 
a stochastic DMC algorithm are set by the step size £ 
and the time increment t, respectively. The key point 
is that for the present application a high spatial and 
temporal resolution is needed only close to the interface, 
whereas far from the dendrite, the coarse-graining length 
and hence the step size can become much larger than 
the fine features of the growing crystal. In practice, we 
choose the step size to be approximately proportional to 
the distance d c b of the walker from the conversion bound- 
ary between the inner (deterministic) and outer (stochas- 
tic) regions, i.e. 



£ w c dcb 



(16) 



with a constant c <C 1. According to Eq. dlif), the 
time increment between updates grows as the square of 
the step size, and hence the walkers far from the den- 
drite have to be updated only rarely. We use dynamical 
lists to efficiently handle the up datin g process, as will be 
described in more detail in Sec. Ill B. For low undercool- 



ings, where the scale of the diffusion field is much larger 
than the dendrite itself and most of the walkers need only 
be updated sporadically, we obtain enormous savings of 
computational time over a straightforward integration of 
the diffusion equation. 

Let us now discuss how the inner and outer regions are 
interfaced. Two essential goals have to be accomplished. 
Firstly, we have to supply a boundary condition at the 
conversion boundary for the integration of the determin- 
istic equations in the inner region, and secondly we need 
to create and absorb walkers at a rate which is propor- 
tional to the local heat flux across this boundary. 

The phase-field equations are integrated in the inner 
region on a regular cubic grid, henceforth called "fine 
grid", with spacing Ax. Each node on this grid contains 
the local values of the phase field ip and the temperature 
field u. We superimpose on this grid another, coarser 
grid, of mesh size L cg = nAx, such that the links of the 
coarse grid intersect the links of the fine grid as shown 
in Fig. |^. The first purpose of this grid is to define the 
geometries of the two simulation regions and of the con- 
version boundary. We describe the "state" of each coarse 
cell by an integer status variable S^g . Here and in the 
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FIG. 2. Sketch of a small part of the conversion boundary 
in two dimensions for n = 4. Each cell of the coarse grid 
(thick lines) contains 16 points of the fine grid (thin lines). 
The fine grid is shown only in the inner region for clarity. The 
shaded cells are conversion cells, and walkers are represented 
by black dots. The boundary F between inner and outer 
regions is indicated by a dashed line. 



following, greek indices (a, (3, 7) label the cells of the 
coarse grid along the x-, y-, and z-directions, whereas 
latin indices (i, j, k) label the nodes of the fine grid. 
All cells which contain at least one node of the fine grid 
where ip > are assigned the status "solid" (S = —2). 
All cells with a center-to-center distance to the nearest 
solid cell smaller than a prescribed length Lb are "buffer 
cells" (S = —1), whereas all other cells belong to the 
outer region. Cells of the outer region which have at 
least one nearest neighbor with buffer status are called 
conversion cells (S — 0) and play the central role in inter- 
facing the two solutions. The dividing surface V between 
inner and outer regions is the union of all the links (or 
plaquettes in three dimensions) of the coarse grid which 
separate conversion from buffer cells (see Fig. ||). Ev- 
idently, as the crystal grows, the geometry of the two 
regions changes, which means that the status variables 
must be periodically updated. Details on this procedure 
are given in Sec. 1IIC. 



We always choose Lb sufficiently large to ensure that 
the phase field is already close to its liquid equilibrium 
value, ip ~ —1, at the conversion boundary. Hence we can 
set tp = — 1 in the entire outer region and treat only the 
standard diffusion equation there. In the initial state, the 
entire system is undercooled to u = — A, and no walkers 
are present. When the crystal grows, it releases latent 
heat which diffuses away from the interface, and hence 
the inner region becomes a heat source for the outer re- 
gion. This heat flux is converted into walkers, each walker 
representing a certain discrete amount of heat. We de- 
fine in each coarse cell an integer variable which 
contains the number of walkers being within this cell at 
time t. For a specific heat which is independent of tem- 
perature, the density of walkers is proportional to the 
difference between the actual and the initial tempera- 
tures, i.e. the temperature in the outer region is related 



5 



to the number of walkers by 



l af3j 



= -A 1 



M 



(17) 



where the constant M fixes the number of walkers in a 
cell that corresponds to the melting temperature u = 0. 

The inner region is completely delimited by conversion 
cells. To fix the boundary condition for the integration 
on the fine grid, it is therefore sufficient to set the field 
u on all nodes of the fine grid in each conversion cell to 
the value specified by Eq. (p7|). The diffusion equation is 
then timestepped in the inner region using the standard 
explicit scheme 
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hjk 



l ijk 
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Hj+lk 
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ij — lk 



(18) 



Note that we have omitted for simplicity the source terms 
due to the phase field, which are zero at the conver- 
sion boundary. Seen on a discrete level, this equation 
can be interpreted as a "pipe flow" equation: the local 
change of u is given by the sum of the "flow" through 
all the discrete links ("pipes"), where, for example, the 
"flow" through a link along x during a timestep is given 
by DAt(ui + ijk — Uijk) / (Ax) 2 . For nodes at the bound- 
ary of the inner region, some links cross the conversion 
boundary V , which means that there is exchange of heat 
with the neighboring conversion cell. This heat flux is 
collected by the conversion cell and stored in a heat reser- 
voir variable H^^. A symbolic manner to describe the 
updating of is 



ff ^7* ~~ H °>0i + (Ax) 2 ( E utgTid U * c 



(19) 



where the sum runs over all the bonds of the fine grid that 
cross r', Ug r id is the temperature on a node of the fine grid 
and u cc is the temperature in the conversion cell given 
by Eq. (refcoarsetemp) . For example, for a conversion 
cell (a, (3, 7) in contact with a buffer cell (a — 1, /3, 7), we 
have (we recall that the linear dimension of a coarse cell 
is L cg — nAx): 



TTt + At _ TTt 

n a0 1 — n ul3~i 



E 



DAt 



^ ^ (Ax) 

j=(f3-l)n+l fc=( 7 -l)n+l v 1 

with i = (a — l)n + 1. 



2 ( U i-ljk U ijk) 



(20) 



If the stored quantity of heat exceeds a critical value H c 
given by 



i d A 
AT' 



(21) 



a walker is created at the center of the conversion cell 
and H c is subtracted from H a p 7 . Conversely, if the lo- 
cal heat flux is negative (heat is locally flowing towards 
the dendrite) and H a p-y falls below — H c , a walker is re- 
moved from the cell and H c is added to the reservoir. 
This algorithm exactly conserves the total heat if the 
contributions of the fine grid, the reservoir variables and 
the walkers are added. In dimensional quantities, each 
walker is equivalent to an amount of heat AQ equal to 



AQ = 



L(nAx) d A 



M 



(22) 



The walkers are restricted to the outer region. If a 
walker attempts to jump across the conversion bound- 
ary the move is discarded and the walker stays at its 
old position until the next update. If c in Eq. ( |l6| ) is 
small enough, such jumps are attempted almost only by 
walkers close to the conversion boundary. Accordingly, 
this procedure is a convenient way of implementing the 
re-absorption of walkers: if a walker stays in a conversion 
cell, the heat flux is more likely to be directed towards the 
inner region, which increases the chances for the walker 
to be absorbed. An alternative method, namely to de- 
posit all the heat contained in a walker in the fine grid 
and remove the walker upon its crossing of the bound- 
ary would create stronger temperature fluctuations on 
the fine grid close to the conversion boundary. 

In summary, the conversion process is handled using 
three auxiliary fields on the coarse grid: the status field 
S a fjy which encodes the geometry of the buffer layer and 
the conversion boundary, the field m a( a 7 that contains the 
number of walkers in each cell and is zero in the inner re- 
gion, and the heat reservoir field iZ Ql a 7 , which is different 
from zero only in conversion cells. Let us comment on 
the size of the grids and the resulting memory usage. The 
fine grid needs to be large enough to accommodate the 
dendrite and the liquid buffer layer during the whole time 
of the simulation. Especially in three dimensions, the re- 
strictions on storage space make it necessary to fully use 
the fine grid. The coarse grid needs to cover at least the 
same space region as the fine grid. As will be detailed 
below, for an efficient handling of the walkers close to the 
conversion boundary, it is desirable to always have some 
portion of coarse grid in front of the conversion boundary, 
and hence the coarse grid should actually cover a slightly 
larger region of space than the fine grid. Since the coarse 
grid has far less nodes than the fine grid (1 node of coarse 
grid for n d nodes of fine grid) , this does not significantly 
increase the storage requirement. In addition, we need 
an array to store the positions of the walkers. The latter 
are represented by "continuous" positions and need no 
grid for their evolution. The walkers can therefore leave 
the region of space where the grids are defined and dif- 
fuse arbitrarily far away from the dendrite, allowing us 
to simulate growth into an infinite medium. The most 
storage-intensive part is the fine grid. In fact, the limit- 
ing factor for most of our three-dimensional simulations 
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is not so much computation time, but rather the storage 
space needed to accommodate large dendrites. 

Finally, let us describe how the different parts of the 
algorithm are connected. The program runs through the 
following steps: 

1. Setup (or update) the status field S a /3-y on the 
coarse grid to fix the geometry of the conversion 
boundary 

2. Calculate the temperature in each conversion cell 
and set the boundary condition for the inner region 
on the fine grid 

3. Timestep the phase-field equations on the fine grid 
and calculate the heat flux between the inner region 
and the conversion cells 

4. Update the heat reservoir variables H a/ 3^ and cre- 
ate or absorb walkers in the conversion cells 

5. Advance the walkers 

6. Repeat steps 2 through 5. From time to time, ex- 
tract the shape of the dendrite and store it for fu- 
ture processing. If the phase boundary has moved 
by more than a coarse cell size, go back to step 1. 

In the following subsections, we will give more details 
on some features of our implementation, such as the up- 
dating of the walkers, the updating of the geometry, the 
choice of parameters, and parallelization. 



B. Updating random walkers 

Before going into details, let us briefly point out similari- 
ties and differences between our method and other DMC 
algorithms. Such methods are widespread in Quantum 
Monte Carlo (QMC) calculations where they are used to 
solve the Schroedinger equation in imaginary time 0. 
Each walker represents a configuration in a usually high- 
dimensional Hilbert space, and the density of walkers is 
proportional to the square amplitude of the wave func- 
tion. In contrast, in our method the walkers evolve in 
real space, and their density represents the temperature 
field. The most important difference, however, is that in 
QMC all walkers are usually updated at the same time, 
whereas in our method some walkers are updated much 
more rarely than others. Therefore, it would be very in- 
efficient to visit every walker in each timestep. Instead, 
we work with dynamical lists. 

To simplify the bookkeeping of the different update 
times, we enforce that updating takes place only at the 
discrete times when the fine grid is updated, i.e. for 
t = iAt, i = 1,2,.... Then, we can make a list for ev- 
ery timestep containing all the walkers that have to be 
updated at that moment. However, these lists greatly 
vary in length and can therefore not easily be accom- 
modated in standard arrays of variables. Therefore, we 
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FIG. 3. Sketch illustrating the configuration of the dynam- 
ical walker lists. Each box stands for a walker, and the full 
arrows indicate pointer variables; the "backbone" array of 
pointers is represented by the downward arrow on the left. 
At time t, walkers are updated and prepended to the lists 
corresponding to their next update time, as indicated by bro- 
ken arrows. 



define a data structure that contains the coordinates of 
one walker plus a pointer variable. Within a given list, 
the pointer associated with one walker indicates the next 
element of the list, or contains an end of list tag if the 
corresponding walker is the last one of the list. An ar- 
ray of pointer variables indicates for each timestep the 
first element of the corresponding list. This array is the 
"backbone" of the list structure. It is easy to add new 
walkers to a list: the pointer of the new walker is set to 
the former first clement of the list, and the pointer of the 
backbone is set to the new walker (see Fig. ||). Lists of 
arbitrary length can be constructed, and every walker is 
visited only when it actually has to be updated. 

At a given time t, the program works through the cor- 
responding list of walkers. The treatment of each walker 
starts by looking up the status of the coarse grid cell cor- 
responding to its position. If the walker is inside a buffer 
cell because the conversion boundary has moved since its 
last update, it is removed. This removal does not violate 
heat conservation because the heat associated with the 
walker is accounted for in the initialization of the tem- 
perature field inside newly created buffer cells (see Sec. 
[IIC, Eq. (24) below). If it is inside a conversion cell, 
and the corresponding reservoir variable < —H c , 

the walker is removed and H c is added to the reservoir. 
In all other cases, the jump distance £ and correspond- 
ing time increment are determined and a new position is 
selected according to Eq. (|l4|). To apply Eq. (|l^) for 
the jump distance i, we need to determine the distance 
of a walker to the conversion boundary. It would be very 
inefficient to calculate this distance for each walker sep- 
arately, especially when the shape of the boundary be- 
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comes complex. Therefore, we use the status field S^g 
on the coarse grid in the outer region to store an approx- 
imate value for this distance, which can then be easily 
looked up by each walker before a jump. Some more 



details are given in Sec. Ill C 



As mentioned above, we restrict the walker updating 
to a discrete set of times. Therefore, the time increment r 
in Eq. (|lj) has to be an integer multiple of the time step 
At, which would not be the case if we directly applied 
Eqs. jltf ) and ([l5]). We solve this problem by defining a 
lower cutoff for the jump distances, 



yj2Dn t At 
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where rit is a fixed integer, and replace the jump dis- 
tances I found from Eq. (|16|) by the closest integer mul- 
tiple of £ m in- We also define a maximum jump length 
l max , mainly to limit the size of the backbone pointer 
array: with a maximum jump distance ^ ma x, each walker 
is at least updated every £^ ax / (2D At) timesteps. Con- 
sequently, the discrete time modulo this number can be 
used to index the pointer variables in the backbone array. 

It should be mentioned that in our list structure, it is 
difficult to find a walker which is close to a given position, 
because all sublists must be searched. This is important 
because the number of walkers in the conversion cells has 
to be known for the interfacing with the inner solution. 
To avoid time-consuming sweeps through the walker lists, 
we update the walker number field m^g on the coarse 
grid whenever a walker jumps. 



C. Updating the geometry 

We now describe more in detail how the status field on 
the coarse grid is setup and adapted to the changing 
geometry. When the dendrite grows, the configuration 
of the buffer layer and the conversion boundary has to 
change in order to maintain a constant thickness Lb of 
the buffer layer. Cells which are part of the outer region 
at the beginning of the simulation may become conver- 
sion cells, then part of the buffer layer, and finally part of 
the dendrite. Under the conditions we want to simulate, 
the crystal may locally melt back, but no large regions of 
space will undergo the transition from solid to liquid, and 
hence we do not consider the inverse status change (from 
buffer to conversion cell, for example). Typically, at low 
undercoolings a readjustment of the geometry becomes 
necessary only after 1000 to 10000 timesteps. Therefore, 
the efficiency requirements are not as stringent as in the 
other parts of the program. 

The procedure starts with a sweep through the fine 
grid. Every cell of the coarse grid which contains at least 
one node of the fine grid where ip > is assigned the 
status "solid" (5*^ = —2). Next, the solid cells at the 
boundary of the dendrite (i.e. each solid cell which has 
at least one neighboring cell which is not solid) are used 



to define the buffer region: all cells with a center-to- 
center distance less than of a boundary cell which are 
not solid are assigned the status "buffer" (S^g = — 1). 
When a conversion cell or a cell of the outer region be- 
comes a buffer cell, we need to define the initial values of 
the two fields on the fine grid. The phase field is set to 
its liquid value, ip = — 1. The temperature is calculated 
from the total heat contained in the cell, taking into ac- 
count both the walkers and the heat reservoir variables 
in the conversion cells in order to ensure that the total 
amount of heat remains conserved, i.e. 



M 



L a/3-f 
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All nodes of the fine grid within the new buffer cell are 
initially assigned this value. The walkers contained in 
the cell are removed. 

All cells of the outer region which are adjacent to the 
buffer, i.e. which have at least one neighbor with buffer 
status, are conversion cells (5 l ^ /37 = 0). When a cell of the 
outer region becomes a conversion cell, its heat reservoir 
variable is initialized at zero. 

Finally, in the outer region, which is comprised of all 
the other cells, the status field is used to store an approx- 
imate value for the distance from the conversion bound- 
ary. A precise determination of this distance is rather 
costly in computation time, because for each cell in the 
outer region, we must calculate the distance to all con- 
version cells and retain the minimum value. A much 
cheaper, albeit approximate method is the following. As 
mentioned, in a conversion cell we have S aS — 0. We 
assign to all cells adjacent to a conversion cell the value 
S t a p 1 = 1. Neighbors of the latter receive the value 
S^p = 2, and we continue this process outward by as- 
signing the value S^p = i + 1 to all cells adjacent to 
a cell with = i. For a relatively simple geometry 

such as a single growing dendrite, the status field can be 
correctly set up on the whole lattice during a single out- 
ward sweep, starting from the center of the dendrite. The 
number assigned to a given cell can be used as a mea- 
sure for the distance. Note that the exact relationship 
of the number to the distance depends on the direction 
with respect to the axes of the coarse grid; our numer- 
ical tests below show, however, that this anisotropy in 
the distance function does not significantly influence the 
dendrite shapes. 

If we follow this procedure, the coarse grid needs to 
cover the entire region of space where the jump distance 
varies. Even though we introduce a large-scale cutoff 
f max , this would become prohibitive in terms of memory 
usage for truly multi-scale problems. Fortunately, such 
a sophisticated scheme for the determination of the dis- 
tance is mainly needed close to the dendrite (for example, 
a walker that enters in the space between two dendrite 
arms needs to make small steps). Once a walker has left 
the vicinity of the dendrite, this rather complicated es- 
timate for the distance to the conversion boundary can 
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be replaced by a simpler one, for example the distance to 
the closest dendrite tip. In consequence, the coarse grid 
needs to cover only a slightly larger region of space than 
the fine grid. 

Finally, let us comment on the integration of the phase- 
field equations in the inner region. We need to know 
which part of the fine grid must be timestepped. This 
information is encoded in the status field S a p on the 
coarse grid. It would, however, be rather inefficient in 
terms of memory access time to integrate the inner re- 
gion "coarse cell by coarse cell". Instead, integration 
proceeds along the spatial direction corresponding to suc- 
cessive memory locations, which is the x-direction in our 
implementation. During the updating of the status field, 
the program determines for each y and z coordinate the 
range(s) to be integrated along x and keeps this informa- 
tion in a lookup table. This table is updated every time 
the status field changes. 



D. Choice of computational parameters 

There are number of parameters in our algorithm which 
can be adjusted to maximize the computational effi- 
ciency. However, certain restrictions apply. Firstly, there 
are various length scales. In order of increasing magni- 
tude, those are: 

1. the lattice spacing of the fine grid, Ax, 

2. the minimum jump length of the walkers, £ m i n , 

3. the size of a coarse-grained cell, L cg = nAx, and 

4. the buffer thickness Lb- 

The minimum jump length should be of the order of the 
inner grid spacing to assure a precise interfacing between 
inner and outer solutions. On the other hand, a larger 
^min means less frequent walker updating. We usually 
worked with £ m [ n « 2 Ax, or n t ~ 10 in Eq. (23). On 



the other hand, 4nin has to be smaller than L cg in order 
to achieve a well-defined coarse-graining. The coarse- 
graining length, in turn, is limited by geometrical con- 
straints. The conversion boundary appears "jagged" on 
the scale of L cg (see Fig. ||). In order to render the ef- 
fects of this coarse geometry irrelevant for the interface 
evolution, the buffer thickness must be much larger than 
this scale, Lb ^S> L cg . We found that L cg w 0.1 Lb is suffi- 
cient to achieve this goal. In our simulations, we mostly 
worked with n = 4 (L cg ss 2^ m ; n ) and n = 8 for larger 
buffer sizes. 

Next, consider the constant of proportionality c be- 
tween the walker jump length and the distance to the 
conversion boundary, d c b- Since the Gaussian random 
vector £ in Eq. (jl4|) has no cutoff, steps of arbitrary 
length are possible, and hence even a walker which is far 
away can jump directly to the conversion boundary. The 



number of such events has to be kept small, because oth- 
erwise the conversion process is influenced by the far field 
with its coarse length and time scales. This goal can be 
naturally achieved by choosing c small enough. For ex- 
ample, for c = 0.1, only jumps with a length of more than 
10 standard deviations can reach the conversion bound- 
ary, which represents a negligible fraction. On the other 
hand, the increase of £ with distance determines the effi- 
ciency of the algorithm, and hence c should be chosen as 
large as possible. We usually worked with c = 0.1, which 
seems to provide a good compromise. 

Finally, the parameter M determines the number of 
walkers per coarse cell and hence the precision of the 
stochastic representation for the temperature field and 
the diffusion equation. Considering Eq. (|l7|), we see 
that the temperature at the boundary of the inner region 
takes discrete values spaced by A/M. In addition, for a 
homogeneous distribution of walkers in a system at u = 0, 
the temperature fluctuations are of order A/\/~M. On 
the other hand, increasing M means longer computation 
time because more random walks have to be performed. 
In addition, the total number of walkers TV necessary to 
simulate a dendrite of final volume V is 



N 



MV 



(nAx) d A' 
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which means that high values of M become prohibitive, 
especially at low undercooling. Fortunately, a good preci- 
sion of the solution can be obtained also by increasing Lb, 
as will be described in Sec. IV. In practice, we worked 



with values of M ranging between 25 and 100. 



E. Boundary conditions and symmetries 

For a two-dimensional dendrite seeded at the origin and 
with arms growing along the x- and y-directions, the sim- 
ulations can be accelerated by taking advantage of the 
cubic symmetry. There are several symmetry axes, and 
consequently it is sufficient to integrate the equations in 
a part of the plane while imposing reflective boundary 
conditions at the proper axes to enforce the symmetry. 
These boundary conditions have to be imposed both on 
the fine grid and for the walkers. For the symmetry axes 
at x = and y = 0, this can be easily achieved by choos- 
ing one of the nodes of the coarse grid to coincide with 
the origin. Then, the two symmetry axes coincide with 
bonds in the coarse grid. On the fine grid, the nodes out- 
side the simulation domain but adjacent to the bound- 
ary are set to the values of their mirror images inside 
the simulation domain after each timestep. Walkers that 
attempt to cross the boundaries are reflected, i.e. in- 
stead of their "true" final position outside the simulation 
domain, its mirror image with respect to the symmetry 
axis is chosen. Another interpretation of this boundary 
condition for the walkers is to imagine that there exists 
an ensemble of "mirror walkers" which are the images of 
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"mirror walker" 



F. Parallelization 
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FIG. 4. Sketch illustrating the implementation of reflecting 
boundary conditions at the symmetry axis x = y. Shown is 
a cell of the coarse grid (solid lines) on the diagonal x — y 
(dashed line). A walker inside the simulation domain (x > 0, 
< y < x) enters the cell. An accompanying "mirror walker" 
(open circles), the image of the walker with respect to the 
symmetry axis, enters the same cell. 

the walkers inside the simulation domain. When a walker 
jumps outside of the simulation domain, its mirror im- 
age jumps inside, and interchanging the walker and its 
mirror, we just obtain a reflection of the walker at the 
boundary as above. 

The latter view is useful when considering the last sym- 
metry axis, the diagonal x = y. While the boundary 
conditions on the fine grid and for the walkers can be 
implemented as before, the conversion process requires 
special attention, because the symmetry axis does not 
coincide with the boundaries of a coarse cell. When a 
walker enters a coarse cell situated on the diagonal, there 
is an additional "mirror walker" entering the same coarse 
cell (see Fig. ||), and hence the number of walkers m^p 
has to be increased by two (or, equivalently, decreased 
by two if a walker leaves the cell). Similarly, walkers 
are created and absorbed in pairs, which means that 
walker creation in such a cell can occur only when the 
heat reservoir exceeds twice the equivalent of one walker. 
In addition, when calculating the heat flux received by 
conversion cells on the diagonal, both the "real" and the 
"mirror" flux has to be taken into account. It is clear 
that this procedure induces an anisotropy in the conver- 
sion process; our tests showed, however, that its effect is 
undetectable for reasonable buffer thickness. 

In three dimensions, the reduction in computational 
resources is even more dramatic. For example, using the 
symmetry planes y = 0, x = y, andx = z, i.e. integrating 
only the domain x > 0, < y < x, z > x, we need only 
integrate 1/48 of the full space, i.e. one eighth of one 
dendrite arm. The planes x = y and x = z can be 
handled as described above, with the exception of cells 
on the diagonal x = y — z. Such cells actually have only 
1/6 of their volume within the simulation domain, and 
for each walker entering a cell, there are 5 mirror walkers 
to be considered. 



Even though our algorithm is very efficient as will be 
shown below, the demands on computation time and 
RAM storage space rapidly increase when the undercool- 
ing is lowered. Therefore, we have developed a parallel 
version of our code for the Cray T3E at the National 
Energy Research Scientific Computing Center (NERSC), 
using the shared memory library SHMEM. 

We are mainly interested in the development of a sin- 
gle primary dendrite branch. Hence, an efficient method 
of parallelization is to divide the simulation domain in 
"slices" normal to the growth direction, and to distribute 
the slices among the processors. In the inner region, 
the integration of the partial differential equations makes 
it necessary to exchange the boundary values between 
neighboring processors after each timestep. This is a 
standard procedure. The more delicate points are the 
handling of the walkers and the updating of the geome- 
try. 

Each processor stores only the parts of the fine grid it 
has to integrate, along with the values of the status field 
in the whole simulation domain. The latter is necessary 
to correctly handle the walkers. For the walkers which 
are far from the dendrite, the average jump distance may 
become much larger than the thickness of a computa- 
tional slice. But if a walker approaches the conversion 
boundary, the conversion process has to be handled by 
the "local" processor which contains the appropriate part 
of the fine grid. Therefore, the walkers need to be redis- 
tributed after their jumps. We have found it sufficient 
to implement "exchange lists" between neighboring pro- 
cessors, i.e. processors which contain adjacent parts of 
fine grid. If a walker jumps to a position outside of the 
local slice, it is stored in one of two lists, correspond- 
ing to "upward" and "downward" motion. After each 
timestep, these lists are exchanged between neighboring 
processors. As most of the walkers make several small 
steps before reaching the conversion boundary, this pro- 
cedure assures the correct redistribution of walkers with 
insignificantly few errors, which arise in the rare case that 
a walker arrives at the conversion boundary after several 
large jumps. 

The only step of the algorithm which needs massive ex- 
change of data between the processors is the updating of 
the geometry: each processor has to determine locally the 
"solid" part of its computation domain, and this infor- 
mation has to be exchanged in order to correctly setup 
the whole status field on each processor. However, as 
mentioned earlier, the geometry is updated only rarely, 
and therefore this part of the algorithm does not repre- 
sent a significant computational burden. We have found 
that the parallel version of our code showed satisfactory 
execution time scaling when the number of processors is 
increased. 
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TABLE I. Computational parameters for the benchmark 
simulations in two dimensions. 



Quantity 


symbol 


Value 


Interface thickness 


Wo 


1 


Anisotropy 




0.025666, 0.00066( 


Effective Anisotropy 


e 


0.025, 0.0 


Relaxation time 


T() 


1 


Kinetic anisotropy 


X 
04 





Grid spacing 


Ax 


0.4 


Timestep 


At 


0.003 


Diffusion coefficient 


D 


10 


Coupling constant 


A 


15.957 


Capillary length 


d 


0.0554 


Kinetic coefficient 


Pa 





Undercooling 


A 


0.3 


Coarse cell size 


n 


4 


Number of walkers per coarse cell 


M 


50 



IV. NUMERICAL TESTS 

The accuracy of the standard phase-field method has 
been assessed in detail by comparison to boundary in- 
tegral results |29). Therefore, to test the stochastic al- 
gorithm it is sufficient to check its results against di- 
rect simulations of the standard deterministic phase-field 
equations. The most critical questions are whether the 
use of the rather coarse lattice for the conversion intro- 
duces spurious anisotropy, and what is the magnitude of 
the temperature fluctuations generated by the stochastic 
treatment of the far field. The main parameters which 
control both of these effects are the thickness of the buffer 
layer and the number M of walkers generated per coarse 
cell. The boundary condition for the inner region is im- 
posed on a coarse geometry with a cutoff scale of nAx, 
and the temperature at the boundary is a stochastic vari- 
able which changes as walkers are created, absorbed, en- 
ter, or leave a conversion box, and which assumes discrete 
values spaced by A/M. When the buffer layer is much 
larger than the size of a coarse cell, 3> nAx, the field 
is "smoothed out" in space and time by the diffusive dy- 
namics. We expect high spatial and temporal frequencies 
to decay rapidly through the buffer layer, and hence the 
evolution of the interface to become smoother as Lb is 
increased. 

We conducted two-dimensional simulations at an inter- 
mediate undercooling, A = 0.3. At this value of A, the 
standard phase-field method can still be used to simulate 
non-trivial length and time scales of dendritic evolution, 
but the length scale of the diffusion field is large enough 
to provide a serious test for the random walker method, 
i.e. the diffusion length is much larger than the thick- 
ness of the buffer layer. Table | shows the computational 
parameters that were used for these tests. Only the first 
quadrant was simulated, with reflecting boundary condi- 
tions at x = and y = 0. Fig. ||(a) shows a comparison 
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FIG. 5. Comparison of standard (deterministic) phase-field 
and random walker method in two dimensions for A = 0.3 and 
€4 — 0.025. (a) Dendrite shapes, represented by the contour 
line 4> — 0, after 200000 iterations, (b) tip velocity versus 
time. 
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FIG. 6. Comparison of "dendrite" shapes without fourfold 
anisotropy after 500000 iterations. 
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FIG. 7. Variance of temperature fluctuations, (m 2 ), as a 
function of the distance from the conversion boundary for 
two values of the walker parameter M. 



of dendrite shapes obtained from the standard phase-field 
and from our algorithm with different buffer sizes. While 
the shapes slightly differ for Lb/ Ax = 20, the curve for 
Lb/ Ax — 40 is almost ^indistinguishable from the deter- 
ministic shape. Fig. ||(b) shows the velocity of the den- 
drite tip along the x-direction, measured over periods of 
500 iterations, versus time. The fluctuations around the 
deterministic value are much larger for Lb / Ax = 20 than 
for Lb /Ax = 40, and for Lb /Ax — 80 (not shown) the 
curve obtained from the stochastic method is very close 
to the deterministic data. For comparison, the diffusion 
length 2D/v at the end of the run is about 400 Ax. 

A particularly sensitive test for the anisotropy of the 
conversion process is the growth of a circular germ with- 
out anisotropy, because such a germ is unstable against 
even smallest perturbations. This can be clearly seen 
from Fig. ^: even though we completely screen the four- 
fold anisotropy created by the lattice (e| = 0), the weak 
next harmonic of the lattice anisotropy, proportional to 
cos 88, destabilizes the circle and leads to the formation of 
bulges in the (21)- and (12)-directions. For Lb/ Ax = 80, 
the stochastic algorithm perfectly reproduces this trend, 
and we can hence conclude that the anisotropy created 
by the coarse structure of the conversion boundary is neg- 
ligibly small. Note that the diffusion field extends to a 
distance of more than 1000 lattice units at the end of this 
run, which means that the larger part of the simulation 
domain is integrated by the stochastic method. 

To quantify the numerical noise, we performed 2-d sim- 
ulations of the simple diffusion equation in a system of 
N x N lattice sites with N = 160. One half of the sys- 
tem (x < 0) was integrated by the stochastic algorithm, 
whereas in the other half (x > 0) we used a standard Eu- 
ler algorithm. The conversion boundary T' hence coin- 
cides with the y-axis, and there is a single column of con- 
version cells along this axis. We used Ax = 1, At = 0.02, 
D = 1, A = 1, and applied no- flux boundary conditions 
at x = ±N/2 and periodic boundary conditions along 
y. The system was initialized at u — everywhere, i.e. 
in the walker region we randomly placed M walkers in 
each coarse cell. When the walkers evolve, fluctuations 



are created in the deterministic region, which plays the 
role of the buffer layer. We recorded u 2 as a function of 
x and averaged over a time which is long compared to 
the diffusive relaxation time of the system, N 2 /D. The 
results for two different choices of the walker parame- 
ter M are shown in Fig. ^. In an infinite homogeneous 
system filled with walkers, the distribution of the num- 
ber of walkers in a given coarse cell is Poissonian, which 
means that the fluctuations in the walker numbers are of 
order \f~M . If this scaling remains valid for the conver- 
sion cells in our hybrid system, we expect (u 2 ) ~ 1/M 
close to the conversion boundary, which is indeed well 
satisfied. As shown in Fig 0, the variance of the tem- 
perature fluctuations rapidly decreases with the distance 
from the conversion layer - by four orders of magnitude 
over the distance of 80 lattice sites. No simple functional 
dependence of (u 2 ) on x is observed. We expect high 
spatial and temporal frequencies to be rapidly damped. 
A theoretical calculation of (u(x) 2 ) seems possible but 
non-trivial because the random variables which are the 
sources of the fluctuations in the deterministic region are 
correlated in space and time by the exchange of walkers 
through the stochastic region and the diffusion of heat 
through the deterministic region. For our present pur- 
pose, we can draw two important conclusions. Firstly, 
for a reasonable thickness of the buffer layer, fluctua- 
tions are damped by several orders of magnitude. The 
residual fluctuations are much smaller than the thermal 
fluctuations represented by Langevin forces that have to 
be introduced in the equations of motion to observe a 
noticeable sidebranching activity Indeed, for suffi- 
ciently large buffer layers we always observe needle crys- 
tals without sidebranches. Secondly, the fluctuations at 
the solid-liquid interface can be reduced both by increas- 
ing the number of walkers and by increasing the thickness 
of the buffer layer, which allows to accurately simulate 
dendritic evolution with a reasonable number of walkers. 

In Table || we compare the run times of our code on 
a DEC Alpha 533 MHz workstation along with the run 
time of the deterministic phase field reference simulation. 
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TABLE II. Execution times of the benchmark simulations 
for various sets of computational parameters. 





M 


L b /Ax 


CPU time (min) 


Deterministic 






1950 




50 


20 


89 




50 


40 


110 




100 


40 


119 



The gain in computational efficiency is obvious. Increas- 
ing the buffer layer from Lb = 20 Aa; to Lb = 40 Ax 
reduces the amplitude of the temperature fluctuations at 
the solid-liquid interface by more than an order of mag- 
nitude, whereas the computation time increases by only 
25%. Comparing the runs with different values of M, we 
see that the walker part of the program accounts only for 
a small part of the total run time. 

From these results, we can conclude that the compu- 
tational effort that has to be invested to simulate a given 
time increment scales approximately as the size of the 
fine grid region, i.e. as the size of the dendrite. This is a 
major advance with respect to the standard phase-field 
implementation on a uniform grid, where the computa- 
tion time scales with the volume enclosing the diffusion 
field. The spatial and temporal scales of dendritic evo- 
lution that can be simulated with our method are hence 
limited by the integration of the phase-field equations on 
the scale of the dendrite. 

All the data shown so far are for two-dimensional sim- 
ulations. We repeated similar tests in three dimensions 
and obtained comparable results for the quality of the so- 
lution and the efficiency of the code. We will not display 
the details of these comparisons here, but rather show an 
example of a three-dimensional simulation under realis- 
tic conditions to demonstrate that our method is capable 
of yielding quantitative results in a regime of parameters 
that was inaccessible up to now. In Fig. g, we show 
snapshot pictures of a three-dimensional dendrite grow- 
ing at an undercooling of A = 0.1 and for a surface ten- 
sion anisotropy e| = 0.025, which is the value measured 
for PVA |2^| . The other computational parameters are 
W = l,Ax= 0.8, e 4 = 0.0284, r = 0.965, 5 = 0.0364, 
ri = 4, M = 50, L b / Ax = 48, D = 24, At = 0.004, and 
A = 39.6 (giving do = 0.0223, /3q = 0), and the simula- 
tion was started from a homogeneously undercooled melt 
with an initial solid germ of radius r — 2Ax centered at 
the origin. During the run, we recorded the velocity v(t) 
and the radius of curvature p(t) of the dendrite tip. The 
latter was calculated using the method described in Ref. 
p9| . With these two quantities, we can calculate the 
time-dependent tip selection parameter 



a*(t) 



2Dd 



(26) 




W)] v{t) 

The results are shown in Fig. |[ In the initial stage 



FIG. 8. Snapshots of a three-dimensional dendrite at 
A = 0.1 after 60000, 120000, 200000, 300000, and 650000 
timesteps (from top left). 



during which the arms emerge from the initial sphere, 
growth is very rapid. Subsequently, the tips slow down 
while the diffusion field builds up around the crystal. At 
the end of the run, the velocity has almost converged to 
a constant value that is in excellent quantitative agree- 
ment with the velocity predicted by the boundary inte- 
gral solution of the sharp-interface steady-state growth 
equations assuming an axisymmetric surface energy and 
tip shape (i.e. the most accurate numerical implementa- 
tion of solvability theory to date pl[). This velocity is 
also in reasonably good agreement with the velocity pre- 
dicted by the linearized solvability theory of Barbieri and 
Langer |33j ] , even though the actual tip radius in both the 
phase-field simulation and the boundary integral calcula- 
tion differ from the tip radius of the paraboloidal shape 
assumed in this theory. A more detailed discussion of this 
point and the entire steady-state tip morphology can be 
found in Ref. 

Remarkably, the selection constant a* becomes almost 
constant long before the velocity and the tip radius have 
reached their steady-state values. This is in good agree- 
ment with the concepts of solvability theory, which stip- 
ulates that the selection of the tip parameters is gov- 
erned by the balance between the anisotropic surface 
tension and the local diffusion field at the tip. To estab- 
lish the correct local balance, diffusion is necessary only 
over a distance of a few tip radii, whereas the buildup of 
the complete diffusion field around an arm requires heat 
transport over the scale of the diffusion length, D/v. Our 
simulation shows that a* indeed becomes essentially con- 
stant soon after the formation of the primary arms. This 
fact can be used to derive scaling laws for the evolution 
of the dendrite arms at low undercooling during the tran- 



sient that leads to steady-state growth [ 
at the end of the simulations, where t 



17[ . Finally, even 
ic dendrite arms 
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FIG. 9. Tip velocity, tip radius and selection parameter 
versus time for the run of Fig. ^. Arrows mark the times 
of the snapshot pictures. Length and time are rescaled by 
do and dg/D, respectively. The steady-state velocity v s s was 
calculated using a boundary integral method |32] . 



are well developed, no sidebranches are visible. We re- 
peated the same simulation for different thickness of the 
buffer layer, and observed no changes in the morphology. 
Tiny ripples can in fact be seen close to the base of the 
dendrite shaft, but the amplitude of these perturbations 
does not depend on the noise strength. We therefore be- 
lieve that this is rather a deterministic instability due to 
the complicated shape of the dendrite base that a begin- 
ning of noise-induced sidebranching. In summary, there 
are at present no indications that the noise created by 
the walkers has a noticeable effect on the morphological 
evolution. 



V. CONCLUSIONS 

We have presented a new computational approach for 
multi-scale pattern formation in solidification. The 
method is efficient, robust, precise, easy to implement 
in both two and three dimensions, and parallelizable. 
Hence, it constitutes a powerful alternative to state of 
the art adaptive meshing and finite element techniques. 
We have illustrated its usefulness by simulating dendritic 
growth of a pure substance from its undercooled melt in 
an infinite geometry. Due to the fact that only a very 
limited amount of "geometry bookkeeping" is required, 
our method can be easily adapted to other experimental 
settings, such as directional solidification. In addition, 
the DMC algorithm is not limited to the present com- 
bination with the phase-field method, but can be used 
in conjunction with any method to solve the interface 
dynamics, as long as the diffusion equation is explicitly 
solved. The adaptation of our method to other diffusion- 
limited free boundary problems is straightforward; prob- 
lems with several diffusion fields can be handled by in- 
troducing multiple species of walkers. 



In view of the results presented here, there is a real- 
istic prospect for direct simulations of solidification mi- 
crostructures for experimentally relevant control parame- 
ters. An especially interesting prospect is to combine our 
method with a recently developed approach to quantita- 
tively incorporate thermal fluctuations fl3ljl in the phase- 
field model. Such an extension should make it possible to 
test noise-induced sidebranching theories fl6| , |34|| in three 
dimensions and for an undercooling range where detailed 
measurements of sidebranching characteristics are avail- 
able |E|26]]3|j3§]. If thermal noise n the liquid region 
outside the buffer layer turns out to be unimportant for 
sidebranching, the straightforward addition of Langevin 
forces as in Ref. |3l| in the finite-difference region (i.e. 
the buffer region plus the solid) should suffice for this ex- 
tension. In contrast, if the noise from this region is im- 
portant, a method to produce the correct level of noise 
in the walker region will need to be developed. Work 
concerning this issue is currently in progress. 

To conclude, let us comment on some possible exten- 
sions and improvements of our method which will be 
necessary to address certain questions. Firstly, we have 
described the method here using an explicit integration 
scheme on the fine grid in the inner region, which enforces 
rather small time steps. We also tested an alternating 
direction Crank-Nicholson scheme in 2-d, which speeds 
up the calculations but makes it necessary to introduce 
corrective terms at the conversion boundary to guarantee 
the local heat conservation. Secondly, for the moment we 
use the stochastic algorithm only at the exterior of the 
dendrite; for other geometries, such as directional solidi- 
fication where the volumes of solid and liquid are compa- 
rable, it might be useful to introduce a second stochastic 
region in the solid. It would also be desirable to combine 
our algorithm with more efficient memory managing tech- 
niques to overcome the limitations due to storage space. 
Finally, a completely open question is whether it is pos- 
sible to combine our stochastic algorithm with a suitable 
method for simulating hydrodynamic equations. This 
would open the way for studies of the influence of convec- 
tion on dendritic evolution at low undercooling, thereby 
extending in a non-trivial way recent studies that have 
been restricted to a relatively high undercooling regime 



ACKNOWLEDGMENTS 

This research was supported by U.S. DOE Grant 
No. DE-FG02-92ER45471 and benefited from computer 
time at the National Energy Research Scientific Com- 
puting Center (NERSC) at Lawrence Berkeley National 
Laboratory and the Northeastern University Advanced 
Scientific Computation Center (NU-ASCC). We thank 
Youngyih Lee for providing the boundary integral results 
and Flavio Fenton for his help with the visualization. Fig. 
p| was created using Advanced Visual Systems' AVS. 



14 



APPENDIX A: PHASE-FIELD METHOD 

We will briefly outline the main features of the phase- 
field method used for our simulations. More details can 
be found in Ref. |§. 

The starting point is the free energy functional, Eq. 
(^), together with the equations of motion for the phase 
field and the temperature field, Eqs. (jl(]) and (pT|). The 
free energy density in Eq. (||) is chosen to be of the form 

/(^-) = 4 + t + a ^( 1 " 2 t + t)- (A1) 

This function has the shape of a double well, with min- 
ima at ip = ±1 corresponding to the solid and the liquid 
phases, respectively. Here, u is the dimensionless temper- 
ature field, A is a dimensionless coupling constant, and 
the term proportional to u on the RHS of Eq. ( Al ) "tilts" 
the double well in order to favor the solid (liquid) mini- 
mum when the temperature is below (above) the melting 
temperature. The coefficient W(n) of the gradient term 
in the free energy (0) determines the thickness of the dif- 
fuse interface, i.e. the scale on which the phase field varies 
rapidly to connect the two equilibrium values. In addi- 
tion, W is related to the surface tension, and exploiting 
its dependence on the orientation of the interface allows 
to recover the anisotropic surface tension of Eq. (H) by 
choosing 



W(n) = Wo 



7o 



(A2) 



The orientation n is defined in terms of the phase field 
by 



(A3) 



Note that this dependence of W on tp has to be taken 
into account in performing the functional derivative, such 
that the explicit form of Eq. (|l0|) becomes 



T {n)d t ip = [ip- Xu(l - tJj 2 )]{1 - ii 2 ) + V ■ [W(n) 2 f<ip] 



m 2 w^ dw{fl) 



d y \V*ip\ W(h 



d z |W| w(h) 



d(d x f/>) 
dW(h) 

dW{h) 



(A4) 



Next, we need to specify the orientation-dependent re- 
laxation time Tin) of the phase-field. In analogy with 
Eqs. (Ia|) and (H) we choose 



r(n) = r (1 - 3c>4) 



1 



4<5 4 



1 - 36 4 



where 64 is the kinetic anisotropy. 

The phase-field equations can be related to the origi- 
nal free boundary problem by the technique of matched 
asymptotic expansions. Details on this procedure can be 
found in Ref. p9fl . As a result, we obtain expressions for 
the capillary length and the kinetic coefficient in terms 
of the phase-field parameters Wo and r(n): 



/J(n) 



a\ T(n) 
1~Wo~ 



a{Wo 



1 - a 2 A 



Win) 2 
Drift) 



(A6) 



(A7) 



(A5) 



where a\ — 0.8839 and a 2 = 0.6267 are numerical con- 
stants fixed by a solvability condition. There is an impor- 
tant difference between this result and earlier matched 
asymptotic expansions of the phase-field equations, due 
to a different choice of the expansion parameter. If the 
coupling constant A is used as the expansion para mete r, 
the first order in A gives only the first term in Eq. (|A7|), 
while the complete expression is the result of an expan- 
sion to first order in the interface Peclet number, which 
is defined as the ratio of the interface thickness and a 
relevant macroscopic scale of the pattern (local radius 
of curvature or dif fusion length). An important conse- 
quence of Eq. (A7) is that the kinetic coefficient and its 
anisotropy can be set to arbitrary values by a suitable 
choice of A and r(n), and in particular we can achieve 
vanishing kin etics (J3(n) — 0). Note that for a r(n) as 
given by Eq. (A5), the kinetic coefficient cannot be made 
to vanish simultaneously in all directions, but for small 
anisotropics choosing 64 = 2e 4 is a sufficiently accurate 
approximation. Furthermore, the ratio do /Wo can be de- 
creased without changing the kinetics by simultaneously 
increasing A and the diffusivity D. This method dra- 
matically increases the computational efficiency of the 
phase-field approach, because the interface width Wo de- 
termines the grid spacing which must be used for an accu- 
rate numerical solution. For a physical system with fixed 
capillary length do, the number of floating point opera- 
tions necessary to simulate dendritic evolution for some 
fixed time interval and system size scales ~ (do/Wo) d+3 
for the choice of phase-field parameters where the inter- 
face kinetics vanish (i.e. Dt/Wq ~ A ~ Wo/do), where 
d is the spatial dimension p9fl . 

We integrate the phase-field equations on a cubic 
grid with spacing Ax, All spatial derivatives are dis- 
cretized using (As) 2 -accurate finite difference formu- 
las, and timestepping is performed by a standard Eu- 
ler algorithm. The use of a regular grid induces small 
anisotropies in the surface tension and the kinetic co- 
efficient. These effects have been precisely quantified 
in Ref. f29[| . Since the grid has the same symmetry as 
the crystal we want to simulate, the presence of the lat- 
tice simply leads to small shifts in the surface tension 
anisotropy and in the kinetic parameters. For example, 
we obtain an effective surface tension anisotropy e| which 
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is sl ightly smaller than the "bare" value 64. Evidently, 



dendrite tip morphology at low undercooling", cond 



the rse of this method restricts the simulation to crystals 
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with symmetry axes aligned to the lattice, but this is not 
a severe limitation in the present study which focuses on 
the growth of single crystals. 
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